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The Wilcox 2006 stress-u; model (also referred to as WilcoxRSM-w2006) has been im- 
plemented in the NASA Langley code CFL3D and used to study a variety of 2-D and 
3-D configurations. It predicted a variety of basic cases reasonably well, including sec- 
ondary flow in a supersonic rectangular duct. One- and two-equation turbulence models 
that employ the Boussinesq constitutive relation were unable to predict this secondary 
flow accurately because it is driven by normal turbulent stress differences. For the NASA 
trapezoidal wing at high angles of attack, the WilcoxRSM-w2006 model predicted lower 
maximum lift than experiment, similar to results of a two-equation model. 


Nomenclature 



drag coefficient 

C l 

lift coefficient 

Cm 

pitching moment coefficient 

C / 

surface skin friction coefficient 

C p 

surface pressure coefficient 

c 

airfoil chord 

D 

duct height 

H 

step height or hill height 

k 

turbulent kinetic energy 

M 

Mach number 

Re 

Reynolds number 

UcL 

centerline velocity 

C r ef 

reference velocity 

U 

velocity component in the x-direction 

u'v' 

turbulent shear stress 

x,y,z 

Cartesian coordinate directions 

Vo 

y at the body 

a 

angle of attack, deg 
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AU difference between upper and lower stream velocities 

S w vorticity thickness 

5 dissipation per unit mass 

uj specific dissipation rate 


Introduction 

Accurately simulating separated flows is one of the most difficult problems in computational fluid dynam- 
ics. For external aerodynamic flows, most researchers use Reynolds- averaged Navier-Stokes (RANS) with 
linear one- or two-equation turbulence models, but predictions for separated flows are not always satisfac- 
tory. The goal of this investigation is to take a first step toward determining if a differential second-moment 
turbulent Reynolds stress model (RSM) can provide improvements over simpler one- and two-equation tur- 
bulence models for such flows. RSMs do not rely on the assumption of a linear relationship between the 
Reynolds stress and strain tensors. They include the effects of streamline curvature and rotation, and thus 
eliminate the need for the special adjustments typically required in one- and two-equation turbulence models. 
They also include normal stress anisotropy near the wall and thus avoid the need for additional modeling 
near stagnation points. Despite the above positive contributions to the understanding of turbulence, most 
investigators shy away from using differential turbulent stress models because of the increased number of 
equations, increased number of sub-models, and lack of numerical robustness. 

Second-moment RSMs consist of five mean- flow conservation equations, six stress equations, and a length 
scale equation. There is no ^-equation because k is determined from the normal stresses. One of the things 
that distinguishes one stress model from another is the length scale equation. A number of approaches have 
been used in the literature to model the turbulent stress equations. They all ultimately lead to a similar 
form. But the length scale equation has been and remains the real problem in turbulence modeling. We 
will continue having problems with turbulence modeling until we can come up with a better length scale 
equation. For the current work, we have selected the existing Wilcox 2006 stress-o; model 1 (also referred to 
as WilcoxRSM-w2006) because it is Galilean invariant, coordinate system independent, and does not require 
the minimum distance to walls. Furthermore, in our experience (and as documented in Wilcox 1 ) o;-based 
models are typically more robust than e-based models, and tend to work better for wall-bounded flows with 
separation. 

The WilcoxRSM-w2006 model 1 was recently implemented in the NASA Langley CFL3D code 2 as part of a 
NASA Research Announcement (NRA) award. The code was thoroughly tested for a variety of 2-D test cases 
provided on the NASA Langley Turbulence Modeling Resource website, a which has complete documentation 
of grids, test data, results, etc. Several examples of these results are provided here. For comparison, results 
using the one-equation Spalart-Allmaras (SA) 3 model and the two-equation Menter’s shear-stress transport 
(SST) 4 model are also included. Additionally, as a 3-D validation test, the WilcoxRSM-w2006 model has 
been applied to flow over the ONERA M6 wing. 

The main purpose of the paper is to use the WilcoxRSM-w2006 model to predict a complex separated flow 
of practical interest: the NASA trapezoidal wing configuration 5 8 at high angles of attack. This configuration 
was the focus of the first High-Lift Prediction Workshop (HiLiftPW-1). Fully turbulent calculations tended 
to under-predict the lift and the magnitude of the pitching moment. Furthermore, agreement between 
computations and experiment tended to deteriorate at higher angles of attack for many of the HiLiftPW-1 
workshop participants. One of the outcomes from additional studies subsequent to the workshop was that 
transition plays an important role for this configuration at the Reynolds number investigated. 9 12 At the 
current time, our implementation of WilcoxRSM-w2006 does not include any transition capability, so in 
this paper fully turbulent results using the WilcoxRSM-w2006 model are compared with results using fully 
turbulent SA and SST. 


Computational Approach 

The CFL3D code is a cell-centered, structured, upwind-biased RANS code. The turbulence equations 
are solved separate from the conservation equations, allowing them to be interchanged easily. This approach 

a http://turbmodels. larc.nasa.gov, accessed 11/12/2013. 
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requires advancing the conservation equations then the turbulence equations in sequence, with the latter 
solved sequentially using approximate factorization. The WilcoxRSM-w2006 model follows this architecture. 

In CFL3D, the convective term is approximated with third-order upwind-biased spatial differencing, and 
the viscous terms are discretized with second-order central differencing. The flux difference-splitting method 
of Roe 13 is employed to obtain fluxes at the cell faces. Advancement in time is accomplished via backward 
Euler, with an implicit approximate factorization method. The turbulence equations are solved with a first- 
order advection scheme. Destruction and diffusion terms in the turbulence models are treated implicitly, and 
production terms are treated explicitly. Additional details regarding the numerical method can be found in 
Krist et al. 2 


Results 

A number of validation cases taken from the NASA Langley Turbulence Modeling Resource website are 
presented first. They include: 2-D mixing, backward facing step, NACA 4412, periodic hill, supersonic 
square duct, and convex curvature duct. In all 2-D cases, grid studies have been independently performed. 
Here, only results are shown for grids that have been deemed to be fine enough for the purposes of RANS 
validation (i.e., further grid refinement made little noticeable difference in results). Predictions for the 3-D 
ONERA-M6 wing are also included in this section. Then, the NASA trapezoidal wing is investigated, where 
the predictions of the WilcoxRSM-w2006 model are compared with the SA and SST models for high angles 
of attack near stall (28, 30, 32, 34, and 37 degrees). 

Validation Results 

Figure 1 compares the turbulent shear stress predictions of WilcoxRSM-w2006 (labeled as RSM in this and all 
other figures throughout the paper) with SA, SST, and experiment for a 2-D mixing layer. This case is from 
Delville et al. 14 In it, two different parallel streams interact, one at 41.54 m/s and the other at 22.40 m/s. 
The Re was 2900 per mm. Figure la shows a sketch of the configuration. Additional details can be found 
on the Turbulence Modeling Resource website referenced above. The results show that WilcoxRSM-w2006 
agreed better with turbulent shear stress measurements from experiment in the near field while SA and SST 
agreed better with experiment at the downstream stations. None of the models captured the progression of 
turbulence levels with downstream distance correctly. WilcoxRSM-w2006 under-predicted the width of the 
shear layer more than the other models at the downstream stations. This is somewhat expected because the 
WilcoxRSM-w2006 model is based in part on calibrations from the Wilcox k-u two-equation model, which 
is known to under-predict spreading rates of free shear layers. 1 

The backward facing step case is shown in Fig. 2. This case is from Driver and Seegmiller. 15 Flow 
conditions were M = 0.128 and Re = 36000 based on step height. A sketch of the configuration is given 
in Fig. 2a (additional details can be found on the Turbulence Modeling Resource website referenced above). 
Results are shown for pressure coefficient, Cp, the skin friction coefficient, C /, and velocity profiles at four 
stations downstream of the step. For this case, it was not possible to obtain a steady state solution for the 
WilcoxRSM-w2006 model. Instead, a time accurate solution was used. The time-averaged results yielded fair 
agreement with experiment. In this case, the SST model predicted the Cf behavior best of the three models 
tested; SA under-predicted Cf and WilcoxRSM-w2006 over-predicted it downstream of reattachment. The 
bubble length of approximately x/H = 6 (as seen in Fig. 2b where skin friction changes from negative to 
positive) was in good agreement with experiment for all three models. 

The NACA 4412 case of Coles and Wadcock 16 was computed at M = 0.09, Re = 1.52 x 10 6 , a = 13.87°. 
Although there may have been non-negligible wall effects in the experiment, 17 this test case was run in 2-D 
with no tunnel walls in order to be directly comparable with other CFD results that have been conducted 
over the years (for example, this is one of the test cases used by Menter 3 4 when he developed the SST model). 
Figures 3 and 4 compare NACA 4412 predictions with experiment. There was a separated region at the 
trailing edge. The SST model appeared to predict the bubble size and extent the best of the three models; 
both WilcoxRSM-w2006 and SA predicted the bubble to be too small. Cp comparisons were similar for the 
three models. Results were good except near the trailing edge, where all models were in error. The SST model 
predicted velocities best compared to experiment. However, the turbulent shear stress was under-predicted 
in magnitude by all three models. 

The next case considered was a separated periodic hill. For this case, the reference solution was the large 
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eddy simulation (LES) of Frohlich et al., 18 which has been shown to be in good agreement with experiment. 
The computations were performed at M = 0.2 and Re = 10595 based on hill height, with periodic inflow and 
outflow boundary conditions at the hill crest and an imposed additional source term driving the mass flow. 
Comparison is made in Figs. 5 and 6. All three models predicted a separation bubble that was significantly 
too long. Similar results were also seen using SST and WilcoxRSM-w2006 for the periodic hill at a lower 
Reynolds number. 19 This particular type of smooth-body separated flow is known to be difficult for RANS 
models to predict. 20 Figure 6a shows velocity profiles at x/H = 2, which is near the base of the back side of 
the hill, near the middle of the separation bubble from the reference LES solution. Here, all three models 
were separated and velocity profiles were qualitatively similar to experiment. However, as shown in Fig. 6b, 
all three models under-predicted the magnitude of the turbulent shear stress in the separated shear layer 
by more than a factor of two. As described in Rumsey, 20 the underprediction of turbulent shear stress 
magnitude results in too little turbulent mixing, leading to significantly longer bubble length prediction. 

Figures 7 and 8 consider the case of supersonic flow in a straight square duct of Davis and Gessner. 21 For 
this case, M = 3.9 and Re = 5.08 x 10 5 based on channel height. The results show that SA and SST cannot 
predict the circulation in the corner regions. This is due to the fact that these models employ the Boussinesq 
eddy viscosity assumption, which cannot differentiate between normal stresses correctly. The normal stress 
differences induce the circulation patterns, so SA and SST predict almost no in-plane velocity components. 
As a consequence of mis-predicting this important corner flow feature, SA and SST were unable to predict 
velocity or skin friction, as seen in Fig. 8. The WilcoxRSM-w2006 model, on the other hand, was in very 
good agreement with experiment. It is worth noting that a quadratic constitutive relation can be used in 
conjunction with SA and SST (instead of the Boussinesq assumption), to improve results for these models 
on this case. 22 

The convex curvature duct case of Smits et al. 23 is shown in Fig. 9. As sketched in Fig. 9a, this flow 
through a channel experienced a rapid 30-degree bend. The Mach number was 0.093 and Reynolds number 
was 2.67 x 10 5 based on channel height. Additional details can be found on the Turbulence Modeling 
Resource website referenced above. The focus was on the convex wall, where curvature is known to cause a 
suppression of turbulence in the boundary layer. All three models showed good agreement with Cp. Both 
SST and the WilcoxRSM-w2006 models produced similar results for Cf and velocity profiles, while SA was in 
worse agreement with experiment. For turbulent stresses, the WilcoxRSM-w2006 model performed the best, 
in particular predicting significantly diminished levels (in magnitude) at the x = 0.030 and 0.183 stations. 
As stated in Wilcox, 1 “Surface curvature... has a significant effect on structural features of the turbulent 
boundary layer... in the absence of ad hoc modifications such effects cannot be accurately predicted with a 
two-equation model... In principal, stress-transport models display none of these shortcomings.” 

Residual convergence behavior is also shown for the convex curvature case, to give the reader a feel for 
the convergence properties of the WilcoxRSM-w2006 model compared to SA and SST. In all three cases, 
mesh sequencing was employed, with the coarser levels run 4000 and 3000 iterations, respectively. As can be 
seen in Fig. 9f, SA and SST exhibited similar convergence rates, while WilcoxRSM-w2006 tended to converge 
poorly during the beginning iterations on each level. Ultimately, WilcoxRSM-w2006 required approximately 
30% more iterations than the other two models to converge to the same level on the finest grid for this case. 
Per iteration, WilcoxRSM-w2006 ran about 45-50% slower than SA and SST. 

The final validation case to be considered was the ONERA M6 wing. 24 This case was computed (fully 
turbulent) at a = 3.06°, M = 0.84, Re c = 11.72 x 10 6 . At these conditions, there was a small region of 
shock-induced separated flow near the wing tip. The grid used was a somewhat coarse 289 x 65 x 49 grid. 
Although the solutions would no doubt benefit from grid refinement in this case (particularly near the shock), 
additional grids were not readily available to us. Note also that this case was run with a symmetry plane at 
the root of the wing; the experiment included a stand-off and splitter plate at this location. As described on 
the NPARC Alliance Validation Archive, 0 this case is considered a “classic CFD validation case,” included 
in numerous validation studies. It is typically computed as we have done here, using a symmetry plane at 
the root. Figure 10 compares the pressure coefficient with experiment. All three turbulence models were 
consistent with each other, and agreed reasonably well with experiment, particularly at and beyond the 65% 
span station. 

b www.grc.nasa.gov /WWW / wind / valid / m6wing/ m6wing,html, accessed 11/12/ 2013. 
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Trapezoidal Wing Results 

The trapezoidal wing is a three-element wing-body configuration described in detail at the HiLiftPW web- 
site. 0 Computations were performed for configuration 1 (slat at 30 degrees, flap at 25 degrees). The Reynolds 
number based on the mean chord (3.30283 ft.) is 4.2 million and the Mach number, M, is 0.2. The grid 
employed was the medium structured multi-zone SX1 grid (no brackets) taken from the HiLiftPW website. 
A schematic of the wing, which was taken from Ref. 12, is shown in Fig. 11. Because we are particularly 
interested in predictions near maximum lift for this case, WilcoxRSM-w2006 was only run at a = 28° and 
above. 

Note that all solutions presented here were started from either free stream or from the output of the 28 
degree case. Converged results were identical in both cases. This is in contrast to some of the participants 
of the HiLiftPW- 1 workshop, who experienced sensitivity to initial conditions for this configuration at high 
angles of attack. * 5 * Figures 12a - 12c compare Cl vs O', Cl vs. C d, and Cm vs. a using SA, SST, and 
WilcoxRSM-w2006. Above a = 28°, both WilcoxRSM-w2006 and SST were similar to each other, producing 
lower lift than SA. The Cl predicted by WilcoxRSM-w2006 did not drop as low as the other models at 
a = 37°, where experiment indicated a large loss of lift. Overall, of the three models, SA provided the best 
comparison with experiment. This is the same conclusion reached at the HiLiftPW-1 workshop, for models 
run fully turbulent. 

Figures 13a - 13c show the pressure coefficient C p for various positions on the flap at an angle of attack 
of 28 degrees. At the upper surface flap-forward station, Fig. 13a shows that none of the models captured 
the peak pressures near the wing tip. This may be in part a result of poor resolution of the wing tip 
vortex. A finer grid in that region should result in improved predictions. 25 Surface pressure predictions 
using WilcoxRSM-w2006 were generally comparable to those of SA and SST, although it yielded slightly 
lower suction pressures on the flap upper surface at 85% span. 

Figures 14a-14d compare velocity profiles vs. normal distance at four different locations. The test results 
were obtained by Hannon et al. 8 in the 14 by 22 ft. NASA Langley wind tunnel using 7-hole probes. Over 
the main wing, there were no measurements in the boundary layer or in the slat wake. The WilcoxRSM- 
w2006 model performed as well as the SA and SST models on the main wing at 15% span and at 83% span. 
Figure 14c shows that WilcoxRSM-w2006 performed well near the flap in the vicinity of the wall, but it did 
not perform as well as SA and SST in capturing the wakes from the upstream elements. Figure 14d shows 
that all models performed poorly over the aft region of the flap. This is the region where surface streamline 
indicated flow separation. 

Overall, the WilcoxRSM-w2006 predictions for the trapezoidal wing at high angles of attack near stall 
were reasonably consistent with those of SST; both models produced lower lift than SA. This was both 
encouraging and discouraging. It was encouraging because it was possible to obtain converged RSM results 
on such a complex configuration; this opens up possibilities for future turbulence model explorations with 
minimal problems from robustness-related concerns. It was discouraging because the RSM results were no 
better than the less expensive SST model for this case. However, bear in mind that the trapezoidal wing 
case is not ideal for this type of validation. Although SA results agreed better with experiment than SST 
or WilcoxRSM-w2006, we reiterate the fact that these computations were performed on the configuration 
without slat or flap brackets, and were run fully turbulent. The brackets were known to have an important 
influence in some regions, 5 and including transition is known to increase predicted lift for this case. 12 

Concluding Remarks 

The main motivation for this research was to determine whether turbulent stress models can be better 
predictors of flow separation than simpler one- or two-equation models. To take steps toward this goal, 
the WilcoxRSM-w2006 model was implemented in the NASA Langley CFL3D code and used to study a 
variety of problems. The model performed reasonably well in all cases considered, but in general it did not 
perform notably better for separated flow cases than SA or SST. In two cases (backward facing step, NACA 
4412) the SST model performed better than SA and WilcoxRSM-w2006; for the ONERA M6 case all three 
models gave similar results; for the shear layer case WilcoxRSM-w2006 was better in the near field and worse 
downstream than the other two models; and for the separated hill case all three models were very poor. In 
two cases (that did not involve separated flow) the WilcoxRSM-w2006 model was notably superior to SA 

c http://hiliftpw. larc.nasa.gov, accessed 11/12/2013. 
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and SST: the supersonic square duct and the convex curvature duct. This makes sense, because the RSM 
accounts for all components of the turbulent stress tensor, and includes the effects of curvature. SA and 
SST must be sensitized through additional modeling fixes in order to account for these effects, which are of 
particular importance for these two cases. 

The trapezoidal wing results showed that the WilcoxRSM-w2006 model results were generally comparable 
to SST in terms of predicted lift near stall. Pressure coefficient and velocity profile comparisons were fairly 
similar for all models. Inadequate grid resolution in specific regions of the flow field may have had some 
influence. Also, brackets and transition have been shown to be an important factor for the trapezoidal 
wing. Transition is not easy to account for when using “fully turbulent” models such as SA, SST, and 
WilcoxRSM-w2006. Future efforts will focus on integrating WilcoxRSM-w2006 with transition prediction 
capability. 

It is important to recognize that both the SST and WilcoxRSM-w2006 models have been “tuned” by 
their creators in order to be able to better predict particular model problems of interest (including some 
separated flows). Menter’s earlier BSL model, which does not include the shear-stress transport correction, 
is well known to have less tendency to separate (it is more diffusive) than SST. It is our belief that for the 
trapezoidal wing, a model with less tendency to separate (such as BSL) will perform better compared to 
experiment near maximum lift when run fully turbulent. However, clearly RANS models are not universal, 
and tuning for particular flows does not always translate into better behavior in general. For example, 
although SST out-performed SA for some of the simple 2-D validation separated flows shown in this paper, 
the opposite occurred for the trapezoidal wing when run fully turbulent. In other words, each model’s 
“tuning” is a compromise in an attempt to achieve reasonable results across a wide range of problems, but 
it is impossible to achieve optimal behavior for all flows. 

The question arises as to whether the current experience with the WilcoxRSM-w2006 model is encouraging 
or discouraging for future efforts with RSMs. Should second-moment RSMs continue to be explored further, 
in an attempt to better understand complex flows, including flow separation? In our opinion, the answer 
is yes. All turbulence models have made assumptions that resulted in their choices for model constants, 
optimized for particular flow fields of interest. The solution to various flow fields depends on the choice 
of these constants, and in some cases the constants could in fact be better represented by non-constant 
functions. Future refinements may improve model prediction of flow separation, and we believe that the 
use of RSMs helps avoid uncertainties associated with additional simplifications that accompany one- and 
two-equation models. 
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(a) Sketch of the configuration 





Figure 1: 2-D mixing shear layer. 
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(a) Sketch of the configuration 


(b) C P 




(d) Velocity profiles 


Figure 2: Backward facing step. 
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(a) Sketch of the configuration 



(b) Experiment (note dark blue areas below airfoil and near top (c) RSM 

of figure contain no data) 



(d) SA (e) SST 

Figure 3: Streamlines for the NACA 4412 airfoil. 
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(c) Stress profiles 

Figure 4: Surface pressure coefficients, velocity profiles, and stress profiles for the NACA 4412 airfoil. 
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Figure 5: Streamlines for the periodic hill. 




(b) Stress profiles 

Figure 6: Velocity and stress profiles for the periodic hill at x/ H = 2. 
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(a) Experiment 



(b) RSM 



(c) SA (d) SST 

Figure 7: Velocity vectors and magnitude for the supersonic square duct at x/D = 50. 



(a) Velocity profiles along diagonal line of duct 


(b) C f profiles along wall 


Figure 8: Velocity profiles and skin friction coefficient for the supersonic square duct at x/D = 50. 
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(a) Sketch of the configuration (b) C p distribution 






Figure 9: Convex curvature duct. 
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(a) 20% span station 


(b) 44% span station 



(c) 65% span station (d) 80% span station 




Figure 10: Pressure coefficient distributions for the ONERA-M6 wing. 
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vel profiles, 15% span 


surface pressure 
(flapfwd) 



^surface pressure 
(50% span) 


surface pressure 
(85% span) 

Figure 11: View of pressure tap locations (on the surface) and 7-hole-probe locations (normal to the surface) 
over the upper surface of the NASA trapezoidal wing. 
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Figure 12: Trapezoidal wing lift curve, drag polar, and pitching moment curve. 
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(a) Flap forward station spanwise. 


(b) Flap at 50% span. 



Figure 13: Trapezoidal wing pressure coefficients for a = 28 degrees. 
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(b) Main at 83% span. 




(d) Aft flap at 83% span. 

Figure 14: Trapezoidal wing velocity profiles for a = 28 degrees. 
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